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A common problem in physics and engineering is the calculation of the minima of 
energy functionals. The theory of Sobolev gradients provides an efficient method for 
seeking the critical points of such a functional. We apply the method to functionals 
describing coarse-grained Ginzburg-Landau models commonly used in pattern formation 
and ordering processes. 
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1. INTRODUCTION 

Many problems in mathematical physics are formulated in terms of finding crit- 
ical points of energy functionals. The recent theory of Sobolev gradients H pro- 
vides a unified point of view on such problems, both in function spaces and in finite 
dimensional approximations to such problems. The aim of this work is to demon- 
strate the use and efficiency of Sobolev gradient techniques in minimising energy 
functionals associated with Ginzburg-Landau models for studying phase transitions 
in alloys and complex fluids. These equations are prototypical for studying pat- 
tern formation or ordering, such as nucleation and spinodal decomposition, that 
are accompanied by instabilities. We illustrate our work with models A and B in 
the Halperin-Hohenberg taxonomy, in which the coarse-grained field or the order 
parameter {OP) is either not conserved (model A) or conserved (model B) [2]. 

A gradient of a functional gives the direction of greatest change per unit change 
in the argument of the functional. Often overlooked is that the direction of a gra- 
dient strongly depends on how the size of arguments of a functional are measured. 
Functionals of interest in physics, particularly energy functionals, commonly in- 
clude derivatives of the arguments. Such arguments have to be considered large 
if some of its derivatives are large. Theoretical considerations of such functionals 
must take this into account but is often overlooked in numerical approximations. 
The theory of Sobolev gradients JlJ is an organized account of how to choose a 
metric for a finite dimensional problem that matches that required for the corre- 
sponding theoretical problem. It is found that a proper matching leads to gradients 
(Sobolev gradients) which are considerably smoother than those normally used 
The result is that the approach to a minimum energy configuration becomes much 
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more efficient. In fact, the improvement in performance using Sobolev gradients 
becomes infinite as mesh size goes to zero. This paper ihustrates this phenomenon 
in some typical problems of interest in phase separation and pattern formation. 
The layout of the paper is as follows: An introduction to Sobolev gradients in Sec- 
tion 2 is followed by a description in Section 3 of Ginzburg-Landau models and 
how Sobolev gradient techniques may be employed. In Section 4 we compare the 
results for minimization using ordinary gradients (functional derivatives) and in an 
appropriate Sobolev space in 1, 2, and 3 dimensions, with different grid spacings 
and with different boundary conditions. If we consider steepest descent as being a 
time evolution from a higher energy state to a lower energy state, then a theoretical 
bound on how large our time step can be is given by the Courant-Freiderichs-Lewy 
(CFL) condition Beyond this limit, the numerical scheme for steepest descent 
will magnify errors in each step. This implies that for the traditional steepest 
descent schemes the step size will have to be decreased as grid spacing becomes 
finer, the dimension of the problem is increased or the order of the derivatives in 
the problem increases. The Sobolev gradient technique avoids these problems 
When we use ordinary gradients we label our results "L2" runs since ordinary gra- 
dients are closely related to attempts at defining a gradient in L2{fl), the space of 
square integrable functions in a region Q. In Section 5 we consider model A', which 
is model A with a constraint, namely, the average value of the OP is conserved. 
This model has a different Sobolev gradient than model A and is an alternative to 
the Cahn-Hilliard equation (model B) when dynamics is not of interest. We com- 
pare results for minimization in L2 using the Cahn-Hilliard approach to model A' 
minimization in the appropriate Sobolev space in 1, 2, and 3 dimensions with dif- 
ferent grid spacings and with different boundary conditions. In Section 6 we extend 
model A' to models of surfactant systems which have higher order derivatives. For 
the models we have studied, the Sobolev gradient technique becomes increasingly 
atttractive as grid spacing is refined, dimension is increased, or the order of the 
derivatives in the problem becomes higher. 



2. THE SOBOLEV GRADIENT APPROACH 



Sobolev gradients essentially provide an organized numerical procedure of de- 
termining prcconditioncrs. An energy functional can be generically written as: 



J(u) 



F{Du), 



(1) 



where is a domain in Euclidean space, u is a member of an appropriate function 
space and Du is a list of (length n, say) consisting of u and all partial derivatives 
of u which are relevant to the problem at hand, is a function on an appropriate 
Euclidean space. For example, consider to be a rectangular region in i?^, F is a 
function on so that 



F{w, r, s) 



w 

T 



w 

T 



for all numbers r, s and D is the transformation Du — (it, Ux, Uy) for all u on U, 
with well-defined partial derivatives. Equation (1) then takes the form 



Jiu) 



u 

T 



u 

Y 



dr, 
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which is one of the functionals we will deal with in this paper. Returning to our 
general considerations of (1), we perform a first variation 

J'{u)h = / F'{Du)Dh. 
Jn 

At this point we depart from custom and do not integrate by parts to obtain the 
Euler-Lagrange equations. Instead, we write 

J'iu)h = I F'{Du)Dh = {Dh, {VF){Du))L,in)-. (2) 

We note that 

{Dh, -D5}L,(n)3 = / {hg + Kg^ + hygy), 
Jn 

the inner product in the Sobolev space ^ [S]. 

By -/j2(il) we mean the Hilbert space of real functions on il in which 

ll/llL(o) -If 

By H^''^{fl) we mean the subspace of L2{fl) consisting of all / so that the norm 

m'm^Hn) = [if + fl + fy) 
Jn 

is defined. 

We introduce a transformation P which is essential to our presentation. Take 
P to be the orthogonal projection of L2(^^)" onto the subspace of all elements of 
the form Du. Such a transformation can be dealt with in a concrete way computa- 
tionally. From (2), 

J'{u)h = {Dh,{VF){Du))L,(n)'^ - {PDh, {VF)iDu))L,in)^ = {Dh, P{VF)iDu))L,in) 

These are legitimate steps since PDh = Dh and orthogonal projections may be 
passed form one side of an inner product to the other. We need one more inner 
product: 

{g,h)s = {Dg,Dh)L^(^nY. 
In terms of this inner product, 

J'{u)h^{Dh,P{VF){Du))L^^^n)^ = {h,(VsJ){u))s, 

where (VsJ)(u) is defined as the first element in the list 

P{yF){Du). 

The function {S/ sJ)(u) is called the Sobolev gradient of J at the element u. To 
make the above calculations useful the projection P must be presented in a suitable 
form and the relevant details are given later. In a number of previous applications 
of the methods (e.g., transonic flow Ginzburg-Landau functionals for supercon- 
ductivity it has been known that Sobolev gradients give vastly superior results 
to those obtained with ordinary gradients. In what follows, slight variations of the 
above will be used, these variations take into account a variety of boundary and 
other external conditions. 
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3. APPLICATION TO GINZBURG-LANDAU MODELS 



Models A and B are defined by the equations ff = and ^ = respec- 
tively, where J is a free energy functional. The static and dynamical properties of 
these models have been extensively studied, primarily in numerical work related to 
coarsening and growth of domains [3 El- The functional J(u) usually has a poly- 
nomial form that depends on the nature of the phase transition as the coefficient 
of the quadratic term changes sign (as a function of temperature, pressure or some 
other thermodynamic variable). The widely used form with terms in and it* is 
associated with a second order or continuous transition, where there is no jump 
discontinuity such as latent heat. 

We seek to minimize the model A free energy functional 



J{u) 



1 



1 



dr 



over some volume subject to certain boundary conditions. The coefficient k deter- 
mines the energy penalty for interfaces. 

In one dimension the problem can be reformulated as minimization of 



J{uq,ui) 



1 



1 



dx 



subject to the constraint that the L2(Xl) functions uq and ui are of the form 



(uO,Ul) = if, fx) 

for some i?^^'^^(ri) function /. We seek a projection operator that maps (uq, ui) in 
L2{^) X L2{^) to the closest point in the subspace consisting of points of the form 
(/, fx). This is given by minimizing the integral 

/ = / [(/ - '^o)' + (/x - u,)^] dx 
over the interval subject to specified constraints. Minimizing / gives the condition 



(1 - dl)f = uo - dxUi. 
A steepest descent scheme in L2{fl) would be of the form 



u ^ u — AVJ(u), 

where A is some scalar and VJ(u) is the variation of J with respect to u subject 
to boundary conditions. We instead perform a steepest descent in the space where 
the gradient is given by the projection we already found: 

V7 7/ N /I ;:,2-,~l f9J{uQ,Ul) dJ{uQ,Ul) 

VJ{uo,ui) = 1 - 5^) M dx 

\ ouq oui 

This is equivalent to changing the norm of candidate functions from 

f^dx 

to 
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Il/ir = I fdx + I fldx 
again subject to appropriate constraints such as boundary conditions. 

4. RESULTS FOR MODEL A 

In this section we report results for model A in one dimension with periodic and 
Dirichlet boundary conditions. The coefBcient k was set to 1.0 for all the numerical 
trials. For periodic boundary conditions, systems of M nodes with spacing h were 
set up with random initial values for the order parameter u such that the avaerage 
value {u) = 0.05 at t = 0. The final minimum energy configuration should have 
u = 1.0 everywhere. The number of iterations, the largest step A that could be 
used, and the CPU time to obtain u > 0.99 everywhere in the system are noted in 
the Tables. The next three entries in the Tables are the number of iterations, step, 
and CPU time required when using the Sobolcv gradient technique. For Dirichlet 
boundary conditions the order parameter u was set to 0.01 everywhere except at the 
ends where u was fixed at zero. The program was terminated when the magnitude 
of the L2 gradient was less than 10^'^ everywhere in the system. 

When minimizing in L2{^) we note that the largest step size that can be used 
for each minimization step decreases as the grid spacing is halved, as is implied by 
the CFL condition. However, steepest des(X!nt using the Sobolev gradient does not 
suffer from this limitation. At each minimization step we first estimate the usual 
L2 gradient, using finite differences to estimate derivatives. Thus, for model A we 
estimate VF = — u — 'V'^u. Using the Sobolev gradient the energy is minimized by 
a step w — > u — A * VsF, where VsF is the Sobolev gradient we want to use. Thus, 
at each minimization step we need to find the Sobolev gradient, given the usual L2 
gradient. This Sobolev gradient satisfies the linear equation (1 — V^)Vs-F = VF. 
This is solved iteratively. The first time we need to calculate the Sobolev gradient 
we do not have a good initial guess, however, in subsequent iterations the Sobolev 
gradient serves as a good initial guess. The Sobolev gradients vary smoothly as the 
minimization progresses and so an iterative procedure is less costly computationally 
than using a direct solver each time. Since the operator (1 — V^) is symmetric, 
positive definite, we use a conjugate gradient solver. Steepest descent and Jacobi 
iteration result in longer run times. 

Results are reported for a single Dec Alpha EV68 CPU. The difference in codes 
for the L2 minimization and the Sobolev space minimization is that in the case 
of the Sobolev space minimization a call to a solver that estimates the Sobolev 
gradient, given the L2 gradient, is made and then the Sobolev gradient is used 
instead of the L2 gradient. 
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One dimensional model A 



Periodic boundary conditions (BCs): 



Nodes M 


si);u-iiig // 


it t'l'ilt KJllS (-^2) 


s1(>i) A f / 1 


CPUs ( L, ) 


literal ions 


sl(>p A 


CPITs 


2iu 


1.0 


1 Q 

io 


U.OZ 


0.00391 


10 


0.6 




0.0195 


2" 


0.5 




1 1 


0.0127 


10 


0.6 




0.0684 


212 


0.25 


173 


0.030 


0.0859 


10 


0.6 




0.325 


213 


0.125 


665 


0.0077 


0.682 


10 


0.6 




1.08 


214 


0.0625 


2674 


0.0019 


5.87 


10 


0.6 




3.07 


2i5 


0.03125 


10514 


0.00048 


51.22 


10 


0.6 




9.75 


Dirichlet BCs: 
















Nodes M 


spacing h 


iterations (L2) 


step A (L2) 


CPUs (L2) 


iterations 


step A 


CPUs 




2lu 


1.0 


38 


0.32 


0.00586 


30 


0.6 


0.0146 




2" 


0.5 


115 


0.11 


0.0244 


33 


0.6 


0.0361 




212 


0.25 


425 


0.030 


0.166 


52 


0.6 




0.159 




213 


0.125 


1660 


0.0077 


1.32 


136 


0.6 




0.906 




214 


0.0625 


6730 


0.0019 


11.64 


370 


0.6 




5.04 




215 


0.03125 


26643 


0.00048 


105.33 


1029 


0.6 




29.69 





For small systems with large spacings the time taken by the solver negates the 
advantage of being able to use a much larger step A when using a Sobolev gradient. 
However, as the system becomes larger and the spacing finer, the Sobolev gradient 
technique is more efficient. 



Two dimensional model A 



Systems now have nodes. 
Periodic BCs: 



M 


h 


iterations (L2) 


A (L2) 


CPUs (L2) 


iterations 


A 


CPUs 


2^ 


1.00 


27 


0.19 


0.005859 


10 


0.6 


0.0107 


26 


0.50 


90 


0.056 


0.0576 


10 


0.6 


0.0693 


2^ 


0.25 


332 


0.015 


0.939 


10 


0.6 


0.709 


2« 


0.125 


985 


0.0038 


14.58 


10 


0.6 


7.52 


29 


0.0625 


3846 


0.00097 


301 


10 


0.6 


77.7 


Dirichlet BCs: 


M 


h 


iterations (L2) 


A (L2) 


CPUs (L2) 


iterations 


A 


CPUs 


2^ 


1.00 


77 


0.19 


0.0127 


36 


0.6 


0.0263 


26 


0.50 


263 


0.056 


0.15 


39 


0.6 


0.181 


2^ 


0.25 


989 


0.015 


2.58 


83 


0.6 


2.46 


28 


0.125 


3909 


0.0038 


53.09 


207 


0.6 


28.38 


29 


0.0625 


15306 


0.00097 


1210.78 


640 


0.6 


387.78 



Again we note that the finer the spacing the less CPU time the Sobolev gradient 
technique uses in comparison to the usual steepest descent. For model A results in 
two dimensions the same step size A can be used for all spacings h when minimizing 
in the appropriate Sobolev space. The step size for minimization in L2 has to 
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decrease as the spacing is refined, we note that it has to decrease much faster in 
two dimensions than in one. 

Three dimensional model A 



Systems now have nodes. 
Periodic BCs: 



M 


h 


iterations (L2) 


A {L2) 


CPUs (L2) 


iterations 


A 


CPUs 


2^ 


1.00 


36 


0.14 


0.303 


8 


0.6 


0.676 




0.50 


124 


0.40 


7.99 


8 


0.6 


7.55 


2^ 


0.25 


494 


0.010 


429.16 


14 


0.6 


91.64 


Dirichlet BCs: 


M 


h 


iterations {L2) 


A {L2) 


CPUs (L2) 


iterations 


A 


CPUs 


2^ 


1.00 


119 


0.14 


0.857 


41 


0.6 


2.32 


26 


0.50 


417 


0.040 


27.57 


55 


0.6 


25.12 


2^ 


0.25 


2115 


0.010 


1395.67 


171 


0.6 


591.31 



The three dimensional models also show that as the spacing becomes finer it is 
advantageous to use the Sobolev gradient technique. We also note from the pre- 
ceding Tables that as the dimension of the problem increases the Sobolev gradient 
technique becomes more efficient. In one dimension Sobolev gradients are more 
efRcient for a spacing h — 0.25, as compared to three dimensions where they are 
more efRcient for spacing h — 0.5. 

5. CONSERVATION CONSTRAINT 

For Model A type systems the order parameter u is not conserved. A Cahn- 
Hilliard [7] or Model B system which would conserve the order parameter is given 

by 

6j{uy 

Su 

Suppose we wish to find the minima of some Model A type functional and we 
require conservation of the order parameter u during the course of the simulation, 
without regard to the actual dynamics. We can use a second projection operator 
to enforce conservation rather than increase the order of our evolution equation by 
two. 

In order that J udu not change, we need to project our gradient onto the sub- 
space of L2{^) functions with integral zero. This is achieved for a function / by 

f/ 

J J y 

We will use the term model A' for model A with this constraint as we do not solve 
for Model B dynamics. The order parameter u is now taken to be a relative concen- 
tration of two fluids A and B with concentrations pA, Pb, such that p = pa + Pb 
and u = PA -PB ^ 

We use the free energy 
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J = J[^{l-u'')-T+^{l + u) log(l/2 + u/2) + ^{l-u) log(l/2 - u/2) + 
f|V«p] (3) 

This free energy contains the entropy of mixing. Phase separation depends on the 
temperature T. When T is greater than the critical temperature Tc — a/2 the two 
phases mix completely. When T is less than Tc there will be domains of positive 
and negative u. The lower T is, the greater can be the possible maximum values 
of |u| at equilibrium. That is, phase separation between fluids A and B is more 
complete at lower T values. 

The model B approach would result in an increase in the order of the derivatives 
of the evolution scheme by two. Imposing conservation through a projection means 
that this can be avoided. As a result, a Sobolev gradient approach for modeling 
systems with conservation constraints is even more efficient in comparison to the 
usual approach. The step size need not be reduced for finer spacings when using 
a Sobolev gradient scheme. Minimization was performed on systems with random 
initial conditions and (u) = 0.05, and a = 2, T = 0.8, k = 1.0 until the magnitude 
of the L2 gradient was less than 10~^ everywhere. By comparing results in two 
and three dimensions we noticei from the Tables that the Sobolev gradient scheme 
is even more efficient in three dimensions than it is in two when compared to the 
traditional approach. 

Two dimensional binary system with periodic BCs: 



M 


h 


iterations {L2) 


A {L2) 


CPUs {L2) 


iterations 


A 


CPUs 


2' 


1.00 


680 000 


0.027 


50.34 


314 


0.95 


0.433 


26 


0.50 


2 516 565 


0.0018 


740 


645 


0.95 


5.26 


2^ 


0.250 


4 420 185 


0.00012 


5187 


1937 


0.95 


98.63 



Three dimensional binary system with periodic BCs: 



M 


h 


iterations 


A 


CPUs 


iterations 


A 


CPUs 


2^ 


1.00 


418 515 


0.012 


6291 


323 


0.95 


33.23 


26 


0.50 


594 233 


0.00086 


68523 


214 


0.95 


418 



These numerical experiments with model A' demonstrate that it is considerably 
more efficient to use a projection to enforce conservation of the order parameter if 
the final equilibrium configuration is all that is important. 

6. SURFACTANT SYSTEMS 

The addition of a surfactant to an oil-water system can be modeled by allowing 
K become negative |H| in the free energy (3). This favors the presence of interfaces 
between the two components of the mixture and thus mimics the action of surfactant 
in allowing the oil and water to "mix" with the formation of bicontinuous domains 
separating the oil and water. We also add a curvature dependent term for a bending 
energy of the form 

to the binary system free energy. By changing 7 one can change the shape of 
domains from circular to oval. The surfactant model enables us to examine how 
the Sobolev gradient approach and the traditional schemes compare when the order 
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of the derivatives increases. The coefficient 7 was set to 1.0 and other parameters 
and initial conditions were as given in Section 5. 



Two dimensional surfactant system with periodic BCs: 



M 


h 


iterations {L2) 


A (£2) 


CPUs {L2) 


iterations 


A 


CPUs 


2^ 


1.00 


4 853 277 


0.0043 


4 696 


43 234 


0.5 


336 


26 


0.50 


27 103 876 


0.000062 


45 250 


4 798 


0.5 


449 


2^ 


0.250 


96 649 780 


0.00000096 


97 327 


5 450 


0.5 


2038 



Three dimensional surfactant system with periodic BCs 

It is clear that a model B minimization with sixth order derivatives will be much 
slower than using model A'. We report results for the Sobolev gradient technique 
only. 



M 


h 


iterations 


A 


CPUs 


2' 


1.00 


30 320 


0.5 


6 636 


26 


0.50 


55 268 


0.5 


630 839 



7. SUMMARY AND CONCLUSIONS 

We have presented minimization schemes for model A Ginzburg-Landau func- 
tionals based on the Sobolev gradient technique ^ |S] . The Sobolev gradient tech- 
nique is computationally more efficient than the usual steepest descent method as 
the spacing of the numerical grid is made finer, the dimension of the problem is in- 
creased, the order of the derivatives in the functional is increased, or a conservation 
constraint is imposed. Our results indicate that Sobolev gradient techniques may 
offer distinct advantages in certain cases, particularly for problems involving func- 
tional that contain spatial gradients such as strain based elasticity problems [5], 
least square formulations of partial differential equations, and electrostatic prob- 
lems that require solving the Poisson-Boltzmann equation. 

An interesting question is whether there exists an optimal metric with respect 
to which the Sobolev gradient works best. It is an open research problem to try to 
find such an optimal metric, even though the optimal one would likely not make 
a large difference computationally in all cases. An example of where there is a 
great difference is in near-singular problems where a weighted Sobolev gradient, 
weighted with the singularity in question, works vastly better ^Hl- The likely fact 
that we cannot yet find an optimal metric may well be responsible for the nonlinear 
dependence of run time on number of grid points noted in this work. 
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